#ifndef __TYPED_MATH__
#define __TYPED_MATH__


#include <math.h>
#include <iostream>

template <typename N>
N do_nothing(N& a) { return (a + a) * (a + 2*a); }

template <typename N>
struct typeLimits {
private:
	static N tol;
	static N under;
	static N e;
public:
	static inline N getTol() { 
		if (tol != 0) return tol; 
		tol = 1; 
		for (int i=0; i<128; i++) {
			volatile N temp = 1 + tol/2;
			do_nothing(temp);
			if (temp == 1) return tol; 
			else tol /= 2; 
		}
		return tol; 
	}
	static inline N getUnderflow() {
		if (under != 0) return under;
		tol = 0;
		N temp = 1 + getTol(); // a number that changes when denomralized
		under = 1;
		for (int i=0; i<0x80000; i++) {
			volatile N temp2 = temp/2;
			do_nothing(temp2);
			if ( temp2*N(2) != temp) break;
			else { under /= 2; temp /= 2; }
		}
		return under;
	}
	static inline N getE() {
		if (e != 0) return e;
		volatile N ans = 1.0;
		N fact = 1.0;
		for (int i=1; i<256 && ans*i*fact + 1 != ans*i*fact; i++) {
			fact *= i;
			ans += 1/fact;
		}
		return (e = ans);
	}
};
template <typename N> N typeLimits<N>::tol = N(0);
template <typename N> N typeLimits<N>::under = N(0);
template <typename N> N typeLimits<N>::e = N(0);

/// TOLERANCE(type) is defined as 2^k such that ((1 + 2^k) != 1) but ((1 + 2^(k-1)) == 1)
template <typename N> inline N TOLERANCE(N a=0) { return typeLimits<N>::getTol(); }
/// UNDERFLOW(type) is defined as the smallest normalized floating point number
template <typename N> inline N MINNORMAL() { return typeLimits<N>::getUnderflow(); }

template <typename N> inline N ABS(N x) { return x<0?-x:x; }

template <typename N> N SQRT(N x) {
	// Newton's iteration
	if (x == 0) return 0;
	if (x < 0) x *= -1;
	N guess = 1;
	N oldguess = 0;
	for (int i=0; i<128 && guess != oldguess; i++) {
		oldguess = guess;
		guess = (guess + x/guess)/2;
	}
	return guess;
}

template <typename N> N LOG(N x) {
	if (x <= 0) return 0;
	N ans = 0;
	// normalize x to between 1/e and e for faster convergence
	while (x > typeLimits<N>::getE()) {
		x /= typeLimits<N>::getE();
		ans += 1;
	}
	while (x*typeLimits<N>::getE() < 1) {
		x *= typeLimits<N>::getE();
		ans -= 1;
	}
	N y = (x-1)/(x+1);
	x = 2*y;
	y *= y;
	int den = 1;
	do {
		ans += x/den;
		den += 2;
		x *= y;
	} while (ans*den + x != ans*den);
	return ans;
}

template <typename N> inline N SQR(N x) { return x*x; }

template <typename N>
N hypot (N x, N y) {
	N m; // min; max in y
	x = ABS(x);	y = ABS(y);
	if (x < y) { m=x; }	else { m=y; y=x; }
	if (!m) return y;
	x = m/y;
	return y*SQRT(1+x*x);
}



template <> inline float TOLERANCE(float a) { return 1.1920928955078125e-07; }
template <> inline float ABS(float x) { return fabsf(x); }
template <> inline float SQRT(float x) { return sqrtf(x); }
template <> inline float LOG(float x) { return logf(x); }

template <> inline double TOLERANCE(double a) { return 2.2204460492503131e-16; }
template <> inline double ABS(double x) { return fabs(x); }
template <> inline double SQRT(double x) { return sqrt(x); }
template <> inline double LOG(double x) { return log(x); }

#ifdef HAVE_LONG_DOUBLE
template <> inline long double ABS(long double x) { return fabsl(x); }
template <> inline long double SQRT(long double x) { return sqrtl(x); }
template <> inline long double LOG(long double x) { return logl(x); }
#endif


// an in-lined template tool. Can compile "generic<CMP<a,b>::MIN>", etc.
template <int a, int b> struct CMP { enum {MIN=a<b?a:b, MAX=a<b?b:a}; };


#endif
